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Abstract 

Shear stress on blood cells and platelets transported In a turbulent flow dictates the fate and biological activity of these 
cells. We present a theoretical link between energy dissipation in turbulent flows to the shear stress that cells experience 
and show that for the case of physiological turbulent blood flow: (a) the Newtonian assumption is valid, (b) turbulent eddies 
are universal for the most complex of blood flow problems, and (c) shear stress distribution on turbulent blood flows is 
possibly universal. Further we resolve a long standing inconsistency in hemolysis between laminar and turbulent flow using 
the theoretical framework. This work demonstrates that energy dissipation as opposed to bulk shear stress in laminar or 
turbulent blood flow dictates local mechanical environment of blood cells and platelets universally. 
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Introduction 

Turbulence is ubiquitous, and is particularly prominent in 
engineered cardiovascular devices as well as pathophysiological 
blood flow. There is strong evidence that turbulence impacts the 
environments of erythrocytes and platelets on the cellular level 
[1,2] in a manner physically distinct from that known in simple 
laminar shear flows such as in a viscometer. While our physical 
understanding of the structure of turbulence and its universal 
properties has received significant attention in the past [3-7] , there 
is no theory that links the statistical properties of turbulence to 
shear stresses physically experienced by cells transported in whole 
blood flow. Shear stress acting on cell membranes is a critical 
mechanical cue that regulates biological activity [7-10] and 
therefore a theory relating turbulence to shear stress environment 
of cells is necessary. 

In turbulent blood flow, the complex spatio-temporal fluctua- 
tions of shear stress leads to hemolysis and platelet activation 
[11,12]. Such phenomena are critical when designing life saving 
devices such as artificial hearts, ventricular assist devices, stents, 
grafts, and heart valves. The phenomena can further impact 
disease such as in the case of an aortic stenosis or atherosclerosis. 
Current models that predict stress experienced by blood cells are 
purely empirical and based on classic experiments [12] that yield 
paradoxical results under laminar and turbulent conditions [13]. 
For a comprehensive review of studies on turbulent blood flow 
related hemolysis and platelet activation, the reader is directed to 
Refs [12,14-17]. As discussed in Ref [12], despite many of these 
studies, there has not been a solid physically justifiable connection 
made between turbulence and the shear stress acting on blood cells 
and platelets. A strong requirement of a physical theory is that it 
should make a link, in a manner independent of laminar and 



turbulent regimes of flow because the pertinent parameter is flow 
at the length scale of individual cells, where it is considered 
laminar. Another aspect that is important to consider is the notion 
of universality of turbulent structures despite the intermittency 
issue [6] . Do the universality properties of turbulence hold in the 
most complex of blood flows considering Newtonian and non- 
Newtonian properties of blood? If so, wiU the distribution of shear 
stress acting on blood cells and platelets be universal? Here the 
term "universal" is used not to imply homogeneous isotropic 
turbulence (HIT), but rather loosely to emphasize the significant 
agreement between the distributions of instantaneous dissipative 
scales in complex in-homogeneous shear flows with the distribu- 
tions observed in HIT as well as the robust presence of inertial 
range scaling[18-21]. 

The goal of this work is to address the above by introducing a 
unified (in the sense of laminar vs. turbulent) physical theory to: (1) 
identify the relevant dynamic properties of flow that link to the 
predicted shear stress experience by cells based on fundamental 
physical arguments. We achieve this for the special case of blood, 
but the underlying physical arguments may hold for any cell 
suspension constrained to the Newtonian regime of its rheology, 
(2) theoretically consider the relevance of non-Newtonian effects 
on the smallest scale turbulent structures for blood based on order 
of magnitude arguments; and (3) test the "universality" of small 
scale structure of turbulent flows in one of the most complex 
turbulent blood flow problems (flow through a bUeaflet mechanical 
heart valve) and experimentally examine the universality of the 
predicted shear stress distribution, thus testing the robustness of 
the theory. Since we lack the capability to experimentally measure 
shear stress on individual cells transported in a turbulent flow, we 
use the new theoretical framework to resolve inconsistency in 
published hemolysis data in laminar and turbulent pipe flow as a 
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Figure 1. Schematic of blood cells in a turbulent "eddy". 

doi:10.1371/journal.pone.0105357.g001 

surrogate validation of the theory. Finally, we note that the 
theoretical framework presented in this work is not limited only to 
blood cells, but applies to any case of suspended cells that are 
sufficiently smaller than the smallest scales of turbulent motion. 

Methods: Theoretical Construction 

Theoretical development is initially constrained to turbulent 
flows where the smallest turbulent eddies are greater than the size 
of the red blood cell (RBC), i.e. the instantaneous minimum size of 
turbulent eddy is always >O(10 \xm). This is an important limit, as 
it is not possible to physically represent turbulent eddies smaller 
than the size of the biggest cells with the single-phase approxima- 
tion of blood. In fact, the single-phase approximation may very 
well be susceptible to errors even when local eddies are within two 
or three times the size of the red-blood cell, i.e. ~25 |J,m. For now, 
let us accept this as a limitation and later discuss repercussions of 
this assumption. Nevertheless, for turbulent eddies >25 |a.m, the 
single-phase continuum representation of blood may be consid- 
ered valid [12]. It is also important to underscore that using a 
continuum single-phase model of blood cannot be equated to the 
assumption of RBCs passively seeded on to a continuous medium. 



On the contrary, the single-phase representation of blood "lumps" 
all effects of the physical reality into the constitutive rheology of 
blood. 

The theoretical construction is next focused on whether small 
scale (dissipative) turbulent structures may be considered Newto- 
nian or non-Newtonian. This is done through order of magnitude 
arguments combined with existing knowledge about blood 
rheology. Subsequently, we introduce basic arguments with 
respect to the nature of turbulent scales of motion in turbulent 
blood flows expected in the cardiovascular system and artificial 
devices. Finally, physical arguments are introduced that link 
turbulence statistics to the distribution of laminar shear stresses 
acting directiy on cells. 

Newtonian vs Non-Newtonian 

It is well known that whole blood significantly deviates from 
Newtonian behavior when local bulk shear rates are below the 
order of 100 s ' [22—24] or when a vessel diameter is less than 
100 |a.m [25]. Now, consider the instantaneous turbulent dissipa- 
tive length scale, ri, defined such that the local instantaneous 
velocity increment (difference), (5m,j, across two points in space 
separated by is such that the local Reynolds number defined as 
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Figure 2. Shear Stress (tp) as a function of eddy size rj plotted for increasing hematocrit based on the equation derived using 
energy balance. 

doi:10.1371/journal.pone.0105357.g002 



l. This instantaiieous length scale is the definition of the 



local turbulent dissipative scale and may be regarded as an 
"instantaneous" Kolmogorov scale [18-20,24,26,27]. Note that 
no assumptions of local isotropy or homogeneous turbulence are 
being made in this framework. For such a small dissipative length 
scale, which in a way is the true measure of the actual size of the 
smallest eddy, the shear rate within such an eddy can be estimated 

as ~ — - ~ Thus the shear rate within such an eddy is set 
' )•/ rj-^ 

purely by v and rj. Now, non-Newtonian behavior of blood dictates 
that V is dependent on y^, which is physiologically 0(1) cSt. It is 

easy to show that the range of y„ ~ — ^ corresponding to eddies of 

sizes \Ofiin<ii<\OQi.im is 0(lO-.s"') <y,| < 0(lO',s^'). Based 
on these order of magnitude estimates, it follows that instanta- 
neous turbulent eddies in the range 10|im<)7< 100|im in any 
cardiovascular blood flow will always be in the Newtonian regime 
considering > O ( 1 00 5 ^ ' ) . In other words, any consequence or 
damage occurring to the cells within these eddies happens while 
the eddy itself behaves as a Newtonian fluid. 

Scales of Motion 

With the notion of blood being Newtonian for turbulent eddies 
10 |J,m and above, let us examine the scales of turbulent motion for 
the smallest scales, without the assumption of isotropy or 
homogeneity. Recall that turbulent flow is characterized by 
complex spatio-temporal fluctuations, M,, superimposed with the 
mean velocity field, C/,. These spatio-temporal fluctuations have 
been phenomenologicaUy represented as turbulent eddies of 



varying length scales. It is important to clarify that physicaUy 
there are no circular "eddies", but the term eddies only refer to the 
Fourier analogy of representing the fluctuating field as an 
summation of kinetic energy over "eddies" ranging from the 
smallest possible scales to the integral length scale. The turbulent 

kinetic energy (per unit mass), k= - UjUi, is transferred across from 

large to the small scales (so called energy cascade) with viscous 
dissipation occurring at the smallest eddies. The Kolmogorov 

length scale, rjj^ = (v'/s) where v is the kinematic viscosity and 
fiis the average local kinetic energy dissipation rate (per unit mass), 
corresponds to the characteristic eddy size based on e. The velocity 
scale of this eddy equals u,^^ = '^/^1k> '''^d therefore the Reynolds 

number, of the Kolmogorov eddy is unity. In the context of 

blood damage, several studies have focused on the relevance of the 
Kolmogorov eddy to predict the capacity of turbulent flow to 
damage blood cells (reviewed in [12]). While r]j^ in theory 
represents the average size of the dissipative eddies (for the case of 
HIT), it does not, in itself, reflect the smallest eddy in a turbulent 
flow in reahstic turbulent flows. In fact, there exist eddies that are a 
fraction of the Kolmogorov eddy, so called sub-Kolmogorov 
eddies, arising from the inherentiy intermittent nature of the 
instantaneous energy dissipation rate field [5,6] . Intermittency is a 
property of all turbulent flows which broadly refers to the 
fluctuations of the fluctuating velocity gradient tensor 

Sij= - (pUi/ dxj-\-duj I dxj) and consequently for R = 2vSjjSjj. With 

energy dissipation rates proportional to the square of velocity 
gradient tensor, the fluctuations in £ are far greater and intense. 
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Figure 3. Schematic of the valve and points of interest upstream and downstream of the valve. 

doi:10.1371/journal.pone.0105357.g003 



with very large departures from its average. Physically these large 
positive fluctuations in energy dissipation occur due to local vortex 
stretching that generate momentarily intense shear regions over 
very small length scales, much smaller than the Koknogorov scale 
itself In other words, when the energy dissipation rate is locally 
above the mean energy dissipation magnitude, the corresponding 
eddy size is sub-Kolmogorov. Thus, the true depiction of micro- 
environment of cells in realistic turbulent flows corresponds to 
their tumultuous experiences in a spectrum of dissipative eddies 
ranging from sub-Kolmogorov scales all the way to the Taylor 
micro-scale that demarks the largest of the dissipative eddies [6] . 
Recent turbulence literature estimate the smallest sub-Kolmo- 
gorov scales to be of roughly Re^^"^ times smaller than the 

Kolmogorov scale [23,26,28]. Here ReL= is the local 

turbulence Reynolds number defined by the integral length scale, 
L, and the local turbulent kinetic energy. To date, studies 
addressing turbulence in blood flow have not taken into 
consideration the issue of intermittency or the sub-Kolmogorov 
fluctuations thereby ignoring the most damaging aspects of 
turbulent flow. In this study, we wiU consider sub-Koknogorov 
fluctuations using direct measurements through the definition of 
the instantaneous turbulent dissipative length scale, )7[18- 
20,24,26,27], while constraining our analysis only to flows where 
the smallest sub-Kolmogorov 10|Um. 

Local Energy Dissipation and Shear Stress Acting on 
Cells. Based on the above picture of scale distributions without 
invoking assumptions of isotropy or homogeneity, and introducing 
the notion of sub-Kolmogorov scales, it is now possible to consider 



a link between macroscopic properties of blood dynamics or eddy 
scales to the local shear stress acting on blood cells. Figure 1 
illustrates a schematic of RBCs in a hypothetical dissipative scale 

V 

eddy of size i-j. Recall that the velocity scale of this eddy is du„ = - 

n 

to yield the condition of locally unit Reynolds number. Given the 
length scale and velocity scale, the instantaneous energy dissipa- 

tion rate corresponding to this smallest eddy is — - = -j. 

Looking at Figure 1 , let us now assume that almost all of this 
energy dissipation physically occurs through viscous straining of 
plasma fluid between cells and that the rate of energy lost in lysing 
or activating cells is negligible in comparison to heat generation 
within plasma. This is a valid assumption if (a) it is known that a 
very small fraction of cells lyse/ activate per unit time, and (b) the 
strain rate of cell membrane deformation is much smaller (by at 
least an order of magnitude) than the strain rate of plasma 
surrounding the cell. The relative difference in the strain rates of 
the cell membrane and the surrounding plasma represents the 
efficiency of energy transfer from the fluid to strain-energy stored 
into the cell membrane. Both (a) and (b) seem valid as hemolysis/ 
activation in most device and physiological flows occur over 
prolonged duration and the magnitude of free hemoglobin 
released is relatively very small compared to the total hemoglobin. 
In particular, the second condition is valid for cases when blood 
cells are lysed over a prolonged exposure to shearing forces as 
opposed to instantaneous lysis which has been observed for whole 
blood shear stresses >450 Pa [29,30]. Complex cardiovascular 
devices such as the bUeaflet mechanical heart valves imposes shear 
stresses in the range < 15 Pa [14]. Furthermore, (b) is supported by 
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Table 1. Basic Turbulence Parameters Upstream. 







y/R 


k (m^/s^) 


L(nm) 




tlo (nm) 


Hk (urn) 


Acceleration 


-0.04 


8.06 E-03 


521 


13.35 


81 


75 




-0.15 


8.81 E-03 


561 


15.05 


80 


73 




-0.25 


7.83 E-03 


580 


14.66 


84 


77 




-0.36 


9.24E-03 


485 


13.33 


75 


70 


Peak 


-0.04 


5.25E-03 


416 


8.61 


88 


83 




-0.15 


476E-03 


283 


5.58 


82 


78 




-0.25 


5.39E-03 


338 


7.09 


82 


78 




-0.36 


5.50E-03 


303 


6.42 


79 


75 


Deceleration 


-0.04 


8.38E-03 


512 


13.38 


79 


73 




-0.15 


9.04E-03 


447 


12.14 


74 


69 




-0.25 


7.40 E-03 


491 


12.07 


81 


76 




-0.36 


7.43 E-03 


595 


14.64 


86 


79 


Regurgitation 


-0.04 


7.53E-03 


365 


9.06 


74 


70 




-0.15 


9.1 1 E-03 


285 


7.77 


65 


61 




-0.25 


7.48 E-03 


346 


8.56 


74 


69 




-0.36 


7.57E-03 


317 


7.89 


72 


67 


Diastole 


-0.04 


9.46 E-02 


7396 


650.03 


70 


57 




-0.15 


9.55 E-02 


4004 


353.50 


59 


49 




-0.25 


9.60 E-02 


7438 


658.48 


70 


57 




-0.36 


5.60E-02 


2746 


185.69 


64 


55 



doi:10.1371/journal.pone.0105357.t001 



the relatively low energy dissipation for both cytoplasm and 
membrane deformation in a RBC tank treading in flow, 
corresponding to 0(10 '-10 ' watts) for shear rates 0(10""), 
which extrapolates to 0(10 watt) for shear rates reaching 
10-' s"' [31,32]. 

Given these assumptions, which appear reasonably valid, the 
rate of energy dissipated within the eddy shown in Figure 1, £, may 
be approximated to the rate of energy dissipated in the fluid that 
surrounds ceUs. The foUowing equation represents the energy 
balance in watts: 



-H)V„ 



Where pg is the density of whole blood, is the volume of the 
eddy, fij, is the dynamic viscosity of plasma, y^is the strain rate in 
the plasma surrounding blood ceUs, and 0<H <l is the fractional 
hematocrit. It is easy to verify that both sides of the equation have 
units of power. In the above equation, if the plasma dynamic 
viscosity is known, then it is straight forward to estimate the viscous 
shear stress, Tp = ^ipYp, within the plasma surrounding the ceUs as a 
fimction of energy dissipation rate. Strictly, the ceUs experience Zp, 

not Tb which can be estimated as %b = HbYb — P'-b^- Clearly, it is 

evident from a physical basis that the relevant dynamical property 
of blood flow that dictates Xp is the local energy dissipation rate, e. 
Recognizing that it is easier to measure jyrather than s, we solve 

equation 1 for Zp = fipfp after substituting s = ^ to give: 



vb 



\-H) 



(2) 



The above equation relates the instantaneous dissipative scale r\ 
in a general turbulent blood flow to the corresponding shear stress 
acting on cells. It is important to note that Zp and Zg are off by a 

1 1/2 



factor Zp/zB- 



Hb{1-H) 



. For typical values, Zp is about 80% 



Local 



of Zb but may easily exceed Zb if locally i/ > 1 — — 

^iB 

variations in hematocrit under turbulent straining may produce 
large discrepancies between Zp and Zb- 

Figure 2 presents the theoretical estimate of shear stress Zp using 
equation 2 over a range of turbxilent eddy sizes Wiim<rj< lO^^w 
and hematocrit 0<i7<0.8. As shown in this figure, Zp ranges 
10 '-10'^ dyne/cm^. We must note here that we extrapolate 
equation 1 for 17 as low as half the size of the RBC. It is interesting 
to note that for an r\ about 5-6 |lm, equation 1 predicts a shear 
stress between 400-500 N/m^, a magnitude well known to 
instantly lyse RBCs [29,30]. This illustrates that the Zp prediction 
of equation 1 up j; approaching 10 |J,m appears to asymptotically 
reach the instantaneous lysis value of 400-500 N/m^ [29,30]. 
Furthermore these shear stresses correspond to the release of 
serotonin from platelets, demonstrating granule release, a process 
involved in platelet activation [33]. 

In real turbulent flows, )j is a dynamically fluctuating quantity in 

space and time, around the statistical measure f/jf = (v^/e) 
Substituting i-jj^ in equation 2 will only represent an average shear 
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Table 2. Basic Turbulence Parameters Downstream. 







y/R 


k (m^/s^) 


L (urn) 




no (fim) 


»lk (uni) 


Acceleration 


0.70 


5.07E-03 


669 


13.61 


102 


94 




0.38 


2.70E-02 


1558 


73.08 


71 


62 




0.27 


2.16E-02 


1543 


64.81 


77 


68 




0.17 


1 .54E-02 


1875 


66.48 


91 


81 




0.06 


9.70E-03 


903 


25.42 


88 


80 




-0.04 


1.23E-02 


1474 


46.70 


93 


83 




-0.15 


2.09E-02 


1651 


68.12 


79 


70 




-0.25 


2.11E-02 


2252 


93.53 


86 


75 




-0.36 


2.16E-02 


1183 


49.68 


71 


63 




-0.67 


1 .60E-02 


484 


17.49 


62 


57 


Peak 


0.70 


3.32E-03 


255 


4.21 


91 


87 




0.38 


1.09E-01 


3537 


333.25 


54 


45 




0.27 


9.38E-02 


2740 


239.76 


53 


45 




0.17 


9.84E-02 


2907 


260.54 


53 


45 




0.06 


5.75E-02 


2212 


151.56 


60 


51 




-0.04 


5.67E-02 


1833 


1 24.73 


57 


49 




-0.15 


1.16E-01 


3873 


376.51 


54 


45 




-0.25 


9.54E-02 


2321 


204.83 


50 


43 




-036 


1.05E-01 


2159 


199.75 


47 


41 




-0.67 


1.13E-02 


423 


12.88 


67 


62 


Deceleration 


0.70 


4.70E-03 


719 


14.10 


107 


99 




0.38 


2.61 E-02 


2009 


92.67 


77 


67 




0.27 


2.30E-02 


2163 


93.81 


82 


72 




0.17 


2.54E-02 


1273 


57.99 


68 


61 




0.06 


2.53E-02 


2004 


91.04 


78 


68 




-0.04 


2.76E-02 


2021 


95.91 


76 


66 




-0.15 


2.33E-02 


1468 


63.99 


74 


65 




-0.25 


2.41 E-02 


1369 


60.69 


71 


63 




-0.36 


2.21 E-02 


1701 


72.27 


78 


69 




-0.67 


1 .70E-02 


1122 


41.76 


76 


68 


Regurgitation 


0.70 


8.42E-03 


269 


7.05 


66 


62 




0.38 


1.00E-02 


10392 


296.89 


172 


145 




0.27 


9.61 E-03 


9143 


256.12 


168 


143 




0.17 


9.17E-03 


12616 


345.11 


188 


158 




0.06 


8.35E-03 


7423 


193.75 


167 


143 




-0.04 


9.18E-03 


4613 


1 26.30 


142 


122 




-0.15 


7.72E-03 


16109 


404.49 


213 


179 




-0.25 


7.63E-03 


8359 


208.61 


179 


152 




-036 


6.40E-03 


6622 


15136 


178 


153 




-0.67 


2.23E-03 


888 


11.98 


148 


138 


Diastole 


0.70 


3.53E-03 


376 


6.39 


99 


94 




038 


1 .OOE-02 


10392 


8.16 


94 


89 




0.27 


9.61 E-03 


9143 


7.40 


93 


88 




0.17 


9.17E-03 


12616 


6.42 


96 


91 




0.06 


8.35E-03 


7423 


5.33 


98 


93 




-0.04 


9.18E-03 


4613 


4.68 


95 


90 




-0.15 


7.72E-03 


16109 


5.66 


103 


98 




-0.25 


7.63E-03 


8359 


5.96 


105 


100 




-036 


6.40E-03 


6622 


4.24 


103 


99 




-0.67 


2.23E-03 


888 


2.15 


110 


108 



doi:1 0.1 371 /journal.pone.01 05357.t002 
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Figure 4. Probability density function Q(r|) during acceleration, peak, and deceleration phases at the points of interest upstream 
(top row) and downstream (middle row). Closure and mid-diastole phases for points of interest upstream are shown in the bottom row. 
doi:10.1371/journal.pone.0105357.g004 



stress acting on RBCs and is perhaps insufficient to reflect the 
intense fluctuations of the dissipation rate that would correspond 
to sub-Kolmogorov scale eddies. Thus, to fially capture the 
statistical nature of shear stress acting on RBC membranes, it is 
necessary to characterize the probability density function of Tp 
denoted P{ip). P(T:p) is related to the probability density 
function of rj, Q{ti), through conservation of probability as: 

P(Tp) = 2()7) — — . Recent turbulence literature indicates that iQ('y) 

may be universal despite the highly intermittent fluctuations of the 
instantaneous dissipation rate field [18,19]. Analytical forms for 
this universal distribution exist with good agreement for experi- 
ments and simulations [24,27]. The key point relevant here is that 
P(i;^)may be universal in the strongly turbulent regions of blood 



flow through devices, and estimated using a simple change of 
variables as: 



P(rp)-- 



2 ^ P 



-3/2, 



0/2 



(3) 



Where C = vb 



[l-HY 

The above arguments and physical construction not only 
provides a mathematical relationship between )7(as a surrogate for 
fi), and Zp, but also yields P{tp). However, how does this 
relationship link turbulence statistics to T^? The most common 
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Figure 5. Normalized probability density function Q(t|/t|o) at the points of interest upstream and downstream at all 
compared to published isotropic result. 

doi:10.1371/journal.pone.0105357.g005 



phases 



and perhaps paradoxical turbulence statistic in the context of 

blood damage has been the Reynolds stress tensor pUjUj. This has 

been extensively discussed in the past (and key issues reviewed in 

[12]) with the note that while it does appear to relate to predicting 

cell damage, there is no consensus with respect to the physical 

relationship between pUiUj and Zp. An alternate model to estimate 

Xp based on cell relative velocities between adjacent eddies has 

been proposed [12] but in our opinion this lacks the physical 

justification for existence of such intense velocity jumps (more 

severe than the smallest sub-Kolmogorov event) over sub-micron 

scales. Nevertheless, based on the phenomenological arguments 

made earlier with the hypothetical eddy, it is easy to see that it is 

the total energy dissipation rate that directly dictates ip. One 

feasible explanation we can offer to explain the robustness of 

Reynolds shear stress in predicting blood damage is that for 

regions in equilibrium, the average rate of energy dissipation, £ is 

related to Reynolds stress given by the equation (proof in [34]): 

dUi ^ dUi. ^ . ^, . 

£ « — UjUj — — , where — — is the mean strain rate tensor, rlugging 
dXj ' dXj ^ 

this approximation in equation 1 and ensemble averaging it, the 

root mean square of Tp is given by: 



^21/2, 



MpPB 



— UiUj 



dUi 
dxj 



1/2 



(4) 



The above discussion is presented for completeness and offers to 
reconcile the rather paradoxical issue (so far) of why Reynolds 
shear stress has been effective in capturing blood damage potential 
of turbulent blood flows [12]. The reconciliation that we offer 
through equation 4 is that it is not the Reynolds stress itself, but the 
product of the Reynolds stress and the local strain rate that 
determines the energy passed on to the cascade and hence the total 
energy dissipation rate, which as we have shown should ultimately 
set the viscous shear stress acting on blood cell membranes. A 
rough order of magnitude verification of equation 4 can be made 



Ns/m^ and mean strain 
■'lAN/rrP'. This strikingly 



based on the contour scale bars in Ref [14] that experimentally 
quantifies the Reynolds shear stresses in a heart valve flow. By 
setting Reynolds shear stress in equation 4~100 N/m^ (i.e. O(IOO) 
N/m^), and substituting f.ip to be ~ 10 

rate to be 0(1000 s"'), we get xfl^ 
agrees with the range of viscous shear stress reported in Ref [14] 
to be 15 N/m^. This agreement demonstrates further that 
Reynolds shear stress when combined with the mean strain rate 
can relate as an order of magnitude estimate to viscous shear stress 
acting on blood cells. The most comprehensive representation of 
shear stress acting on ceUs however is undoubtedly the hypoth- 
esized universal distribution function P{ip). 

A Test of Universality and Quantification ot P(Xp). With 
the reasonableness of the Newtonian assumption for physiological 
turbulent eddies, we interrogate the most complex turbulent blood 
flow problem — the pulsatile turbulent field surrounding a 
bUeaflet mechanical heart valve using high-resolution phase-locked 
particle image velocimetry to calculate P(Xp) and assess its 
universality. Turbulence in mechanical heart valve flows has 
received enormous attention in medical research with unresolved 
issues in relation to the apphcability of Reynolds shear stresses on 
blood cells [12]. Briefly, a bUeaflet mechanical heart valve was 
experimentally subjected to physiological aortic conditions and the 
instantaneous turbulent velocity field was captured in the vicinity 
of the valve during representative phases of the cardiac cycle. 
Points of interest within the measurement region where complexity 
is expected (Figure 3) were further interrogated to quantify the 
probability density function Q{r]) and the corresponding proba- 
bility density function of the shear stress,P(Tp) for a hematocrit of 
48%. 

Particle Image Velocimetry (PIV) of Bileaflet Mechanical 
Heart Valve. Experiments are similar to that described in Ref 
[35]. The fluid was seeded with polyamide tracer particles (Dantec 
Dynamics, Denmark) distributed in the 5 to 30 microns range. A 
laser sheet illuminated the central plane normal to the b-datum 
line. The laser was generated using the Photonics Industries 
DM40-527 diode-pump Qjswitched laser (Photonics, Bohemia, 
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Figure 6. Probability density function of instantaneous sKiear stress (tp) during acceleration, peak, and deceleration phases at the 
points of interest upstream (top row] and downstream (middle row). Closure and mid-diastole phases for points of interest upstream are 
shown in the bottom row. 
doi:10.1371/journal.pone.0105357.g006 



NY) with optics to covert the output beam into an expanded laser 
sheet. The laser had an initial thickness of approximately 1 mm, 
which was focused down to less than 200 microns within the 
measurement region using a spherical lens (f = 1 m). The valve 
was oriented such that the measurement plane bisected both 
leaflets at the central plane of valve model. A Photron Fastcam 
SA3 CCD high speed video camera (Photron, San Diego, CA) 
synchronized to the laser system via a high speed controller (HSC) 
(La Vision, Ypsilanti, MI) captured focused images of the 
illuminated polyamide particles within the laser sheet in the 
measurement plane. The image area of interest was 1 .5D wide and 
ID tall with the valve body centered. Image distortion due to 



curvature of the acrylic tube was compensated in-situ with a 
calibration plate consisting markers placed in a regular square grid 
with 1 mm spacing. The DaVis calibration algorithm (LaVision 
Inc, Germany) automatically tracks the markers and a map to 
evaluate the corrected image. Corrected image generated of the 
calibration plate verified successful calibration and distortion 
correction. The PIV setup achieved a raw data spatial resolution of 
roughly 27 |J.m/pixel. PIV measurements were performed in 
double-frame mode with a laser pulse separation time. 
At = .500 us. This ensured adequate particle displacements in the 
range of 10-15 pixels thus maximizing the accuracy of instanta- 
neous velocity measurements to within 2% error. 
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Figure 7. Probability density function of instantaneous normalized shear stress (tp/tK) for all points of interest at all phases. Notice 
the close data collapse for Tp/tK above 0.5. The probability of Tp>0.5tk; is about 0.36 (shaded region). 
doi:1 0.1 371/journal.pone.01 05357.g007 



Images were pre-conditioned by first subtracting the minimum 
image from every image acquired. Instantaneous two-dimensional 
velocity field was calculated from the raw particle images using 
cross-correlation processing with a multi-pass scheme. The initial 
interrogation window size for die multi-pass scheme was at 32 x32 
pixels, which progressively reduced to 8x8 pixels. Interrogation 
window overlap was frxed at 50%. Post-processing of the vector 
data included a median filter that rejected vectors outside 3 
standard deviations of the neighbor vector. Gaussian smoothing 
was used to reduce noise in the vector data. An in-house Matlab 
code was developed to post-process these raw velocity measure- 
ments to derive statistical properties, specifically PDFs. 

Peak locking index was calculated to be between 0.02-0.19. 

Peak locking index is defined as 4|^— -| where <j> is the first 

moment of the probability distribution function of the absolute 
fractional distance in pixels to the nearest integer pixel displace- 
ment. If the probability distribution is uniform in the pixel 
displacement range 0 to 0.5, then the pixel locking index is zero, 
indicating no pixel locking. A value of 1 indicates 100% pixel 
locking (i.e. no sub-pixel displacements recorded). Our range of 
0.02-0.19 for peak locking index is far below 0.25, which is the 
threshold for minor peak locking artifacts. 

Error Considerations. The sources of error in our mea- 
surements are due to resolution as well as random error. Random 
errors are addressed in this study through statistical averaging of 
repeated measurements (N = 500) and statistical comparisons. This 
section briefly outlines the errors in accuracy due to limited 
resolution of the measurement techniques at hand. The conser- 
vative error estimate in velocity is <2% (i.e. particle displacements 



may be off by ±0.2 pixels out of total displacement of 10 pixels). 
Laser pulse timing errors are negligible in comparison. 

Validation. In order to check if the valve chamber and the 
acrylic model valve provide equivalent results to clinical quality St. 
Jude Medical valve. Figure 5 of Ref [35] compares noii- 
dimensionalized leaflet kinematics, and the downstream velocity 
profile at x/D = 0.33 during peak systole to published results for a 
clinical quality St. Jude Medical valve [36] . 

Calculation of Q{ri) and P(ip). PDFs of the local dissipation 
scale, Q{r]) was calculated for each interrogation point in Figure 3. 
The approach is similar to that previously described [20,37] with 

the qualifying condition for t] being 0.9 < < 2. Briefly, the 

instantaneous r\ is estimated from the local velocity increment 
between PIV velocity measurement points. We have shown in Ref 
[37] that a significant portion of the local dissipative scale PDF, 
Q{t]) may be derived if the velocity measurement resolves a 
sizeable portion of the dissipative range; i.e. PIV resolution well 
resolves the Taylor microscale. Given that the PDF is derived from 
the histogram of the occurrence of (y, and that the measured 
variable in the above inequality is 5u^, it is straight forward to 
propagate the percent error in the instantaneous velocity of the 
PIV measurements to the uncertainty in the PDF. Our uncertainty 
of 2% in velocity translates to an uncertainty of 2.8% in &u,f. 
Given that the inequality is the only qualifying criteria, the 
uncertainty in the PDF may be achieved by perturbing the upper 
and lower limits of the inequality. To be conservative, we 
recalculated the PDFs by incorporating a 10% variation in the 
limits and found that the resulting PDFs with this additional 10% 
uncertainty in c)M,j insignificantly influenced the shape of the PDF. 
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Energy Dissipation Rate {Watts) 

Figure 8. A plot of plasma hemoglobin data from Kameneva et. al (2009) with respect to total energy dissipation rate calculated 
from reported wail shear stress and Reynolds numbers. Figure illustrates how total energy dissipation consistently captures blood damage 
irrespective of laminar or turbulent flow regimes. 
doi:10.1371/journal.pone.0105357.g008 

The PDF P(ip) was determined by directly calculating the 
occurrences oitp for every occurrence of (7 using equation 2. Non- 
dimensionalization parameters were calculated for each interro- 
gation point and are listed in Tables 1 and 2. 

Results: Dimensional Pdfs of x\ and Xpfor Flow Near 
a Blleaflet Mechanical Heart Valve 

Figure 4 presents the probability density function Q{r\) at points 
of interest (defined in Figure 3) during specific phases of the 
cardiac cycle. Upstream of the valve, the peak of Q{t]) 
(representing the mode of the fluctuating rf) is consistentiy observed 
to be between 100-200 |J,m during forward flow. The smallest 
eddies captured in the distribution function are in the range 50- 
70 |J,m. The same characteristics are observed upstream of the 
valve during closure phase. However, during mid-diastole, there is 
a significant change in Q{v]) with respect to the varying location of 
points. In particular, except for the furthest location from the 
centerline, the locations closer to the centerline correspond to a 
significant leftward shift of Q{r])- For these locations, the mode of 
r] is in the range 80-90 \xxa, with the smallest r] about 40 |J.m. 
Downstream of the valve, there is significant variation in Q{v]) 
characteristics as a function of the lateral position relative to the 
centerline during acceleration, peak, and deceleration phases 
(Figure 4). The smallest r] is between 40-50 |j,m and the mode oir] 
between 90-200 |im. Positions located within the side orifice jets 
above and below the centerline do not show significantiy different 
Q{yi) characteristics compared to the upstream forward flow 
characteristics. 

The above shifts in Q{vi) characteristics clearly point to the 
shght leftward shift whenever increased turbulence is expected. For 
instance, locations slightly off from the centerline upstream of the 
valve win experience the high shear regions of the regurgitation jet 
during mid-diastole. Similarly, the leftward shifting of Q{y]), also 
coincides with the shear layer locations downstream of the two 
leaflets consistent with the results in a more classical flow problem 



[36]. For all the locations, it appears that the smallest eddies are in 
the neighborhood of about 40 |J,m, which is consistent with the 
literature [12]. On the contrary, these small eddies appear to be 
relatively rare events with most of the dissipative eddies in the 
turbulent zones around 80 |a.m. This may be lower in locations 
where significantly greater turbulence exists. 

Figure 5 presents the normalized probability density function 
Q(l/lo) fo'" each of the positions of interest upstream and 
downstream at different phases of the cardiac cycle. vIq is defined 
as r}o = LReJ^^'^^ where L is the local integral length scale, and 
Re]^ = k^/-L/\< as discussed in Refs. [18,26] is a scale very close to 

rij^ = LRej^^^* . Also presented is Q{>l/rio) ^or the case of 
homogenous isotropic turbulence (HIT) from highly resolved 
direct numerical simulations [26] for comparison. Figure 5 shows 
very good agreement of the experimentally derived PDFs with the 
HIT expectation. The scatter around the HIT expectation may be 
attributed to experimental uncertainty as well as weak dependence 
of Q{r]lr\o) on local mean shearing [20]. Specifically, we have 
shown that this weak dependence occurs when the Corrsin length 

scale, = (s/S'^^) where S is the mean shear magnitude, 

approaches the mean shear-viscosity length scale, (v/S)'^^ [20]. 
Nevertheless, from a practical standpoint these results confirm the 
largely valid assumption of universal small scale structures despite 
the highly pulsatile nature of the turbulent flow past a complex 
device. This is particularly significant given the classical assump- 
tions behind universality of turbulence are highly restrictive (i.e. 
very high Reynolds number, fuUy developed, equilibrium, 
stationary etc.). 

Figure 6 presents the raw PDFs of Ty,. The peak (mode) of the 
distribution for Zp ranges between 5 and 20 dynes/cm^ with the 
right tail extending to about 60 dynes/cm'. Tables 3 and 4 list the 
maximum observed Zp corresponding to the minimum r] at each 
position for all recorded phases for the cardiac cycle. To assess the 
universality of P(Zp), we examine the normalized probability 
density function P(Zp/zp ) (see Figure 7) for the different points of 
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interest and cardiac phases. Xp^ is defined by substituting 

r\j^ = LRe^ '^'^ in equation 2, thus defining a characteristic 
Kolmogorov scale shear stress acting on blood cells as: 

'PK = ]lY^^ReJ (5) 

The normalized PDFs show strong data collapse indicating a 
practical universality which is expected based on the collapse 
observed for Q{ri/rjo)- Note that the maximum observable Zp was 
roughly regardless of location of measurement or cardiac 
phase for the specific problem. The tale of Pizpjtp^) drops 
exponentially around this magnitude. A theoretical limit to the 
maximum shear stress Tp^^^ may be estimated utilizing the 
smallest possible sub-Kolmogorov scale ris = LReL~^ [26,28] to 

Discussion 

Resolving Historically Inconsistent Hemolysis in Laminar 
and Turbulent Pipe Flow 

The argument introduced to relate turbulence properties to the 
shear stress environment is a simple balance of energy dissipation. 
Thus, regardless of whether a flow is laminar or turbulent, 
isotropic or not, energy dissipated must highly correlate to the fate 
of suspended cells. In fact, energy dissipation was classically 
recognized as the hemodynamic parameter that dictates hemolysis 
based on simple physical arguments [38,39]. We test if this can 
resolve one of the paradoxical results published in Ref [1 3] where 
hemolysis was measured in fuUy developed laminar and turbulent 
pipe flows while maintaining the same wall shear stress. Given that 
the total shear stress of the mean flows for both these cases are 
identical, the conventional thinking expects the same levels of 
hemolysis. However, turbulent cases produced significandy higher 
hemolysis [13]. To resolve the lack of an appropriate physical 
interpretation, we present a revised graph (Figure 8) that plots the 
measured plasma free hemoglobin (a measure hemolysis) as a 
function of total energy dissipation rate (in watts) in the pipe based 
on available information of waU shear stress and Reynolds 
numbers. Briefly, based on the wall shear stress magnitudes, 
length of the pipe, and its diameter the pressure drop can be 
calculated using a simple force balance. We also calculated flow 
rates from the given Reynolds number, \'isc{)sit)', and pipe 
diameter information. The pressure drop and flow rate were 
multiplied to yield the total energy dissipation rate in watts. 
Figure 8 confirms that indeed, regardless of whether the flow is 
laminar or turbulent, the total energy dissipated produces a single 
monotonic functional dependence for hemolysis. The functional 
form with the best correlation was an exponential fit with 
R^ = 0.91 compared to linear, and power fit models. Nonetheless, 
this result supports our fundamental assumptions and confirms 
that total energy dissipation is the fundamental parameter that 
dictates the mechanical environment of suspended cells. 

One must note however that in the above pipe flow data, 
exposure time was not a relevant quantity as in all the experiments 
of ref [13], the total exposure time was fixed. Thus the resulting 
trend shown in Figure 8 reflects the exponential growth in 
hemolysis due to the complex instantaneous energ\' dissipation 
"histories" of the turbulent cases. Further studies are necessary to 
develop appropriate blood damage models that take into account 



instantaneous energy dissipation rate histories and exposure times 
for pulsatile apphcations. 

In the Context of a Cellular Environment 

Cells are highly responsive to changes in their energetic 
environment, as demonstrated by the effectiveness of m vivo laser 
injury models [40]. Thermal energy is capable of inducing 
alterations to RBC morphology and can result in hemolysis in 
extreme cases [41,42]. Exacerbating the problem, a cascade of 
events can occur upon RBC lysis including platelet activation 
stimulated from adenosine diphosphate [43]. More direct studies 
of mechanical stress have demonstrated numerous effects on 
suspended blood cells, which can be combined with an energy 
balance equation, as presented here to determine the degree of cell 
activity. High shear stress is known to cause hemolysis and can 
further result in shear-induced platelet activation and shear- 
induced platelet aggregation, which has become known as SIPA, a 
feature that occurs independent of activation agonists. This 
process depends on the magnitude and the duration of shear 
stress [33]. More complex gradients in shear stress have been 
shown to promote activation-independent platelet aggregation 
involving the binding protein, von Willebrand factor (VWF) [44], 
a feature that may be important for the intermittent nature of 
turbulence. Platelets combined with VWF can result in aggrega- 
tion while in suspension, independent of an adhesive surface 
substrate [45] . For sufficient shear stress, VWF can change from a 
globular confirmation to an extended chain [46,47], which can 
result in s(-lf-association, creating a structure that is \-cry adherent 
to platelets and can create a physical barrier for other cells [48]. 
Shear stress in bulk flow or at a surface will also influence the 
inflammatory response. Activated platelets, alone, release a 
number of prothrombotic and proinflammatory molecules 
through a-granules and microparticles [49-.'il]. These processes 
may interact with mechanotransductive mechanisms in mono- 
cytes, which are highly responsive to membrane tension through 
cellular activation, signaling, and migration [52,53]. 

The current model may better predict the local shear stress 
environment for these processes and for dictating the function of 
binding proteins, platelets, RBCs, and monocytes, especially for 
the complex flow environment that can exist in medical devices 
and atherosclerosis [54]. As it becomes clear how cells respond to 
their mechanical environment, it will be increasingly important to 
develop unifying theories, such as the one presented in this paper 
to better determine the mechanical conditions for cells in a 
physiological or pathophysiological environment. It will further be 
critical for medical device development to consider the relation- 
ship between the cellular response and mechanics. For example, it 
has already been shown that self-association of VWF can result in 
acquired von Willebrand Disease in devices, including ventricular 
assist devices [55,56], a blood disorder that can be amplified by 
shear-induced receptor shedding in platelets [57]. 

In addition to predicting cellular functional responses to the 
local mechanical environment, it is also important to consider the 
influence on transport in the fluid flow. The complex biorheology 
of blood is highly dependent on the local shear environment [58] . 
The local environment influences the number of collisions between 
cells and the rotational nature of the cells, both of which are 
largely depcnrk^nt on the hematocrit. Transport in the flow may be 
critical for tliroml)us growth rates and the movement of 
microparticles or cellular agonists [59,60]. 

Limitations 

The proposed theory assumes that the distribution of eddies are 
not largely influenced by factors, such as the mean shear in the 
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flow environment, an assumption sliown to be reasonably valid for 
the mean shear at the points examined. However, there may be 
much stronger shearing in medical devices such as LVADs, where 
strong mean shear effects need to be accounted for to revise the 
distribution function of the dissipative scales of equation 3. In our 
probabilit)' density functions, we did see slight shifts, which may be 
caused by these influences, discussed in detail elsewhere [20,37]. 
However, the shifts were relatively minor in the current study. 

Conclusion 

We introduced a theory to predict the mechanical environment 
of cells in turbulent l)l<)od flows through physical arguments 
applicable to any in-homogenous turbulent flow. It is argued that 
for the case of physiological blood flow, the dissipative turbulent 
eddies are well in the Newtonian regime of blood rheology. We 
experimentally showed that the "universal" prediction of dissipa- 
tive eddy distributions is valid even in a highly complex flow 
generated past a bileaflet mechanical heart valve, and that 
consequendy the shear stress distributions experienced by blood 
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